﻿using System;
using SoftwareConsulting.BI11.NumericalMethods.DhbFunctionEvaluation;
using SoftwareConsulting.BI11.NumericalMethods.Interfaces;
using SoftwareConsulting.BI11.NumericalMethods.Iterations;
using SoftwareConsulting.BI11.Studio.InternalServices;

namespace SoftwareConsulting.BI11.NumericalMethods.Statistics
{
    /// <summary>
    /// Класс гамма-распределения
    /// </summary>
    public class GammaDistribution : ProbabilityDensityFunction
    {	    
        private double _alphaShape;//параметр формы распределения	    
        private double _betaScale;//параметр масштаба распределения	    
        private double _norm;//норма распределения (кэшируется для ускорения)	    
        private double _a;//константа для генератора случайных чисел (кэшируется для ускорения)
        private double _b;//константа для генератора случайных чисел (кэшируется для ускорения)
        private double _q;//константа для генератора случайных чисел (кэшируется для ускорения)
        private double _d;//константа для генератора случайных чисел (кэшируется для ускорения)	    
        private IncompleteGammaFunction _incompleteGammaFunction;//незавершенная гамма-функция распределения (кэшируется для ускорения)

        /// <summary>
        /// Возвращает/устанавливает параметр формы распределения
        /// </summary>
        protected double AlphaShape
        {
            get
            {
                return _alphaShape;
            }
            set
            {
                _alphaShape = value;
            }
        }

        /// <summary>
        /// Конструктор класса по умолчанию
        /// </summary>
        protected GammaDistribution()
        {
        }

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="shape">форма распределения</param>
        /// <param name="scale">масштаб распределения</param>
        public GammaDistribution(double shape, double scale)
        {
            if (shape <= 0)
                throw new ArgumentOutOfRangeException("shape");
            if (scale <= 0)
                throw new ArgumentOutOfRangeException("scale");
            SetParameters(shape, scale);
        }

        /// <summary>
        /// Конструктор класса
        /// </summary>
        /// <param name="histogram">Гистограмма, из которой вытягиваются параметры по наилучшим предположениям</param>
        public GammaDistribution(IHistogram histogram)
        {
            if (histogram.Minimum < 0)
                throw new ArgumentOutOfRangeException();
            double shape = histogram.Average;
            if (shape <= 0)
                throw new ArgumentOutOfRangeException();
            double scale = histogram.Variance / shape;
            if (scale <= 0)
                throw new ArgumentOutOfRangeException();
            SetParameters(shape / scale, scale);
        }

        /// <summary>
        /// Возвращает среднее распределения
        /// </summary>        
        public override double Average
        {
            get 
            { 
                return _alphaShape * _betaScale; 
            }
            set 
            {
                Assert.Fail(Assert.Id.Assert00282);                
            }
        }

        /// <summary>
        /// Задает новые параметры (не проверяя их на корректность)
        /// </summary>
        /// <param name="shape">форма распределения</param>
        /// <param name="scale">масштаб распределения</param>
        public void SetParameters(double shape, double scale)
        {
            _alphaShape = shape;
            _betaScale = scale;
            _norm = Math.Log(_betaScale) * _alphaShape + GammaFunction.LogGamma(_alphaShape);
            if (_alphaShape < 1)
                _b = (Math.E + _alphaShape) / Math.E;
            else 
                if (_alphaShape > 1)
                {
                    _a = Math.Sqrt(2 * _alphaShape - 1);
                    _b = _alphaShape - Math.Log(4.0);
                    _q = _alphaShape + 1 / _a;
                    _d = 1 + Math.Log(4.5);
                }
            _incompleteGammaFunction = null;
        }

        /// <summary>
        /// Возвращает значение функции плотности распределения
        /// </summary>
        /// <param name="x"></param>
        /// <returns>вероятность, нахождения случайной переменной меньше или равной, чем X</returns>
        public override double GetDistributionValue(double x)
        {
            return IncompleteGammaFunction.GetValue(x / _betaScale);
        }

        /// <summary>
        /// Возвращает незавершенную гамма-функцию
        /// </summary>
        private IncompleteGammaFunction IncompleteGammaFunction
        {
            get
            {
                if (_incompleteGammaFunction == null)
                    _incompleteGammaFunction = new IncompleteGammaFunction(_alphaShape);
                return _incompleteGammaFunction;
            }
        }

        /// <summary>
        /// Возвращает эксцесс распределения
        /// </summary>
        public override double Kurtosis
        {
            get 
            { 
                return 6 / _alphaShape; 
            }
        }

        /// <summary>
        /// Возвращает название распределения
        /// </summary>
        public override string Name
        {
            get 
            { 
                return "Gamma distribution"; 
            }
        }

        /// <summary>
        /// Возвращает/устанавливает параметры распределения
        /// </summary>
        public override double[] Parameters
        {
            get 
            { 
                return new double[] { _alphaShape, _betaScale }; 
            }
            set 
            { 
                SetParameters(value[0], value[1]); 
            }
        }

        /// <summary>
        /// Возвращает случайное вещественное число в соответствие с распределением
        /// </summary>
        /// <returns></returns>
        public override double GetRandomValue()
        {
            double r = _alphaShape > 1
                        ? GetRandomValueForAlphaShapeGreaterThan1()
                        : (_alphaShape < 1 ? GetRandomValueForAlphaShapeLessThan1() : GetRandomValueForAlphaShapeEqual1());
            return r * _betaScale;
        }

        /// <summary>
        /// Возвращает случайное вещественное число для формы равной 1
        /// </summary>
        /// <returns></returns>
        private double GetRandomValueForAlphaShapeEqual1()
        {
            return -Math.Log(1 - RandomGenerator.NextDouble());
        }

        /// <summary>
        /// Возвращает случайное вещественное число для формы больше 1
        /// </summary>
        /// <returns></returns>
        private double GetRandomValueForAlphaShapeGreaterThan1()
        {
            int maxTries = short.MaxValue;
            for (int i = 0; i < maxTries; i++)
            {
                double u1 = RandomGenerator.NextDouble();
                double u2 = RandomGenerator.NextDouble();
                double v = _a * Math.Log(u1 / (1 - u1));
                double y = _alphaShape * Math.Exp(v);
                double z = u1 * u1 * u2;
                double w = _b + _q * v - y;
                if (w + _d - 4.5 * z >= 0 || w >= Math.Log(z))
                    return y;
            }
            Assert.Fail(Assert.Id.Assert00283);
            return 0;
        }

        /// <summary>
        /// Возвращает случайное вещественное число для формы меньше 1
        /// </summary>
        /// <returns></returns>
        private double GetRandomValueForAlphaShapeLessThan1()
        {
            double y = 0;
            int maxTries = short.MaxValue;
            for (int i = 0; i < maxTries; i++)
            {
                double p = RandomGenerator.NextDouble() * _b;
                if (p > 1)
                {
                    y = -Math.Log((_b - p) / _alphaShape);
                    if (RandomGenerator.NextDouble() <= Math.Pow(y, _alphaShape - 1))
                        return y;
                }
                y = Math.Pow(p, 1 / _alphaShape);
                if (RandomGenerator.NextDouble() <= Math.Exp(-y))
                    return y;
            }
            Assert.Fail(Assert.Id.Assert00284);
            return 0;
        }

        /// <summary>
        /// Возвращает эксцесс распределения
        /// </summary>
        public override double Skewness
        {
            get 
            { 
                return 2 / Math.Sqrt(_alphaShape); 
            }
        }

        /// <summary>
        /// Возвращает как строку
        /// </summary>
        /// <returns></returns>
        public override string ToString()
        {
            return string.Format("{0} ({1:####0.00000},{2:####0.00000}", Name, _alphaShape, _betaScale);
        }

        /// <summary>
        /// Вычисляет значение, основываясь на функции плотности
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        public override double GetValue(double x)
        {
            return x > 0 ? Math.Exp(Math.Log(x) * (_alphaShape - 1) - x / _betaScale - _norm) : 0;
        }

        /// <summary>
        /// Возвращает дисперсию распределения
        /// </summary>
        public override double Variance
        {
            get 
            { 
                return _alphaShape * _betaScale * _betaScale; 
            }
        }
    }
}
